#include <vtkSmartPointer.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>

#include <vtkLookupTable.h>
#include <vtkFloatArray.h>
#include <vtkCellData.h>
#include <vtkPolyData.h>
#include <vtkPlaneSource.h>
#include <vtkSphereSource.h>

#include <vtkNamedColors.h>
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkImageData.h>
#include <vtkSphereSource.h>
#include <vtkMetaImageWriter.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkImageStencil.h>
#include <vtkPointData.h>
#include <vtkCutter.h>
#include <vtkPlane.h>
#include <vtkStripper.h>
#include <vtkLinearExtrusionFilter.h>
#include <vtkSTLReader.h>

#include <vtkFeatureEdges.h>
#include <vtkProperty.h>
#include <vtkClipPolyData.h>

#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkImageData.h>
#include <vtkSphereSource.h>
#include <vtkMetaImageWriter.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkImageStencil.h>
#include <vtkBMPWriter.h>
#include <vtkPointData.h>
#include <vtkCubeSource.h>

#include <QFileInfo>

#include <algorithm>
/// 将文件转到 mdh -> dicom

/// 计算两点间的距离
static double distance(const double a[3], const double b[3]) {
    return sqrt((b[0]-a[0])*(b[0]-a[0]) + (b[1]-a[1])*(b[1]-a[1]) + (b[2]-a[2])*(b[2]-a[2]));
}

bool writeToMhdFile(const char* stlPathName, char *resultMhdFileName)
{
    QFileInfo fi = QFileInfo(stlPathName);
    if (!fi.exists()) {
        return false;
    }
    printf("fileName:%s, baseName:%s, suffix:%s\n",
           fi.fileName().toStdString().c_str(),
           fi.baseName().toStdString().c_str(),
           fi.suffix().toStdString().c_str());

    vtkSmartPointer<vtkSTLReader> reader = vtkSmartPointer<vtkSTLReader>::New();
    reader->SetFileName(stlPathName);
    reader->Update();

    //vtkSmartPointer<vtkCubeSource> cubeSource = vtkSmartPointer<vtkCubeSource>::New();
    vtkSmartPointer<vtkPolyData> pd = reader->GetOutput();
    vtkSmartPointer<vtkImageData> whiteImage = vtkSmartPointer<vtkImageData>::New();

    double bounds[6];
    pd->GetBounds(bounds);
    bounds[0] -= 2.0;
    bounds[1] += 2.0;
    bounds[2] -= 2.0;
    bounds[3] += 2.0;
    bounds[4] -= 2.0;
    bounds[5] += 2.0;

    double spacing[3]; // desired volume spacing
    spacing[0] = 0.2;
    spacing[1] = 0.2;
    spacing[2] = 0.2;
    whiteImage->SetSpacing(spacing);

    // compute dimensions
    int dim[3];
    for (int i = 0; i < 3; i++) {
        dim[i] = static_cast<int>(ceil((bounds[i * 2 + 1] - bounds[i * 2]) / spacing[i]));
    }
    whiteImage->SetDimensions(dim);
    whiteImage->SetExtent(0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1);

    double origin[3];
    origin[0] = bounds[0];
    origin[1] = bounds[2];
    origin[2] = bounds[4];
    // 这样和原始图偏移了 （-2，-2，-2）
    whiteImage->SetOrigin(origin);
    whiteImage->AllocateScalars(VTK_UNSIGNED_SHORT,1);

    // fill the image with foreground voxels:
    unsigned char inval = 255;
    unsigned char outval = 0;
    vtkIdType count = whiteImage->GetNumberOfPoints();
    for (vtkIdType i = 0; i < count; ++i) {
        whiteImage->GetPointData()->GetScalars()->SetTuple1(i, inval);
    }

    // 写个球, 无内容的写不出来
    // short *pData = (short *)whiteImage->GetScalarPointer();


    std::vector<std::vector<double>> ironBalls = {
        {17.5, 60, 29},
        {17.5, 120, 29},
        //{83.5, 120, 29},
        //{102.5, 60, 29},
    };

    //    for(uint n = 0; n < ironBalls.size(); n++) {
    //        double x = ironBalls[n][0]+2.0;
    //        double y = ironBalls[n][1];
    //        double z = ironBalls[n][2];

    //        for(int k = z/spacing[2] - 20; k < z/spacing[2] + 20; k++) {
    //            for(int j = y/spacing[1] - 20; j < y/spacing[1] + 20; j++) {
    //                for(int i = x/spacing[0] - 20; i < x/spacing[0] + 20; i++) {
    //                    double p0[3] = {i*spacing[0], j*spacing[1], k*spacing[2]};
    //                    if(distance(p0, ironBalls[n].data()) < 1.1) {
    //                        whiteImage->GetPointData()->GetScalars()->SetTuple1(k * dim[0] * dim[1] + j * dim[0] + i, 10000);
    //                    }
    //                }
    //            }
    //        }
    //    }

    // polygonal data --> image stencil:
    vtkSmartPointer<vtkPolyDataToImageStencil> pol2stenc =
      vtkSmartPointer<vtkPolyDataToImageStencil>::New();
    pol2stenc->SetInputData(pd);
    pol2stenc->SetOutputOrigin(origin);
    pol2stenc->SetOutputSpacing(spacing);
    pol2stenc->SetOutputWholeExtent(whiteImage->GetExtent());
    pol2stenc->Update();

    // cut the corresponding white image and set the background:
    vtkSmartPointer<vtkImageStencil> imgstenc =
      vtkSmartPointer<vtkImageStencil>::New();
    imgstenc->SetInputData(whiteImage);
    imgstenc->SetStencilConnection(pol2stenc->GetOutputPort());
    imgstenc->ReverseStencilOff();
    imgstenc->SetBackgroundValue(outval);
    imgstenc->Update();

    vtkSmartPointer<vtkMetaImageWriter> writer =
      vtkSmartPointer<vtkMetaImageWriter>::New();
    QString mhdName = fi.baseName() + ".mhd";
    sprintf(resultMhdFileName, "%s", mhdName.toStdString().c_str());
    writer->SetFileName(resultMhdFileName);
    writer->SetInputData(imgstenc->GetOutput());
    writer->Write();

    return true;
}

void showUsage(char *name)
{
    printf("Usage: %s STL_file DICOM_folder\n", name);
    printf("       requirement:\n");
    printf("       compile InsightToolkit and copy ImageReadDicomSeriesWrite to /usr/bin or /usr/local/bin\n");
}

int main(int argc, char** argv)
{
    if (argc != 3) {
        showUsage(argv[0]);
        return -1;
    }

    char resultMhdFileName[128];
    bool ret = writeToMhdFile(argv[1], resultMhdFileName);
    if (!ret) {
        showUsage(argv[0]);
        return -2;
    }

    char cmdLine[256];
    sprintf(cmdLine, "ImageReadDicomSeriesWrite %s %s", resultMhdFileName, argv[2]);
    system(cmdLine);
}
